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Abstract 

Thermodynamic and magnetic properties of Ising models denned on the triangular Apollonian 
network are investigated. This and other similar networks are inspired by the problem of covering 
an Euclidian domain with circles of maximal radii. Maps for the thermodynamic functions in 
two subsequent generations of the construction of the network are obtained by formulating the 
problem in terms of transfer matrices. Numerical iteration of this set of maps leads to exact 
values for the thermodynamic properties of the model. Different choices for the coupling constants 
between only nearest neighbors along the lattice are taken into account. For both ferromagnetic 
and anti-ferromagnetic constants, long range magnetic ordering is obtained. With exception of 
a size dependent effective critical behavior of the correlation length, no evidence of asymptotic 
criticality was detected. 
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I. INTRODUCTION 



The investigation of magnetic models on scale invariant networks has attracted the atten- 
tion of scientists since the . Besides the fact that, on such graphs, renormalization 
procedures can lead to exact results^, they have been explored as models for systems that 
are not translational invariant, neither in the positions of the spins nor in the coupling con- 
stants mediating the interactions between them. In this respect, the analysis of disordered 
and aperiodic models on scale invariant graphs, which include hierarchical lattices 0, 0, 0], 
Cayley trees [?J or Sierpinski gaskets and carpets^, Q], have provided valuable insight into 
the behavior of critical phenomena of non homogeneous systems on Euclidean lattices. 

A further family of scale invariant graphs are Apollonian networks, the simplest of which 
is illustrated in Fig. 1. This lattice can be defined based of the ancient problem of filling 
space with spheres, first tackled by the Greek mathematician Apollonius of Perga |l0|. In 
its two dimensional version, corresponding to the problem of the plane filled by circles, the 
nodes of this network are defined by the positions of the centers of the circles, while edges are 
drawn between any pair of nodes corresponding to pairs of touchin g ci rclesjlll] . The resulting 
network, corresponds to the contact force network of the packing[12]. Apollonian networks 
can also be used to describe generically other scale-free situations like space-filling porous 
media or the connections between densely located cities on which one is interested in 
fluid flow, car traffic or electric supply. Therefore it is useful to study not only the geometric 
properties of Apollonian networks but also transport and ordering on them. 

In a previous paper ^j], we analyzed several of its properties. In particular, we have 
shown that it has small-world properties, a scale-free degree distribution, a very high clus- 
tering coefficient and a very short diameter, all this having been confirmed independently 
by Doye |l5[. Moreover, the Apollonian network can be embedded in the Euclidian plane, 
what is not the case for other scale invariant lattices, e.g. hierarchical lattices or Ca yley 
trees. One can of course define other similar lattices based on modified packing rules |lfi| . 
Although the numerical values for exponents depend on the topology of each realization, 
basic properties characterizing complex networks remain the same. 

In our first work, we have also devoted our attention to several physical models on the 
Apollonian network (electrical resistance, percolation, magnetic ordering), pointing out the 
most striking features. In this work, we review our investigations on the properties of several 
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Ising models on the Apollonian network, and present a thorough discussion, regarding on one 
hand some of the details of the Transfer Matrix (TM) methods and, on the other hand, the 
most important thermodynamic and magnetic properties. Several choices for rules defining 
the values of the coupling constants are considered, including both ferro (F)- and anti- 
ferromagnetic (AF) interactions. As will be shown, we found long range magnetic ordering 
for almost all choices of couplings without any noticeable evidence of a phase transition to 
a paramagnetic phase at a finite temperature. For particular choices of F and AF bonds, 
the geometry of the network induces the presence of competition and frustration within 
closed loops of odd number of sites, giving rise to residual entropy and changes in the 
correlation length. The behavior of this quantity deserves a detailed discussion, as it points 
to a transition from long to short range correlation only for finite size systems, in a similar 
way as discussed for magnetic models on scale-free lattices 17 1 

The rest of this paper is organized as follows: in Section II we discuss our models. 
In Section III we obtain the maps for the free energy and its derivatives, as well as for 
the correlation length. Results are discussed in Section IV, while concluding remarks are 
presented in Section V. 



II. ISING MODELS 



The Apollonian network is constructed recursively. In each generation, it incorporates 
a new set of sites, which correspond to the centers of the new circles added to the packing 
filling the holes left in the previous generation. In the present work we consider the lattice 
which starts with three touching circles drawn on the vertices of an equilateral triangle, and 
the packing problem is restricted to filling the space bounded by these three initial circles, as 
shown in Figure la. If n denotes the current generation of the network, the number of sites 
N(n) is asymptotically three times that of the previous generation, i.e. N(n+1) = 3N(n) — 5, 
or N(n) = (3 n_1 + 5)/2. The number of edges linking nodes increases with n according to 
B(n+1) = B(n)+3(N(n+l)-N(n)). As a consequence, B(n) = (3+3")/2, B(n)/N(n) -> 3 
in the limit of large n, so that on average, each site is linked to six other sites, which is the 
coordination number of the triangular lattice. 

Once the Apollonian network has been defined, it is possible to define many different 
models on it. In this work we focus on a set of interacting Ising spins = ±1 placed on 
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FIG. 1: Third generation of construction of the Apollonian network. Sites represented by squares, 
circles and diamonds are introduced in the first, second and third generation. Links represented 
by dotted, dashed and solid lines correspond to n = 1,2 and 3 in eq|HJ). 

each site of the network. Interactions are restricted to pairs of spins placed on nodes linked 
by edges, as described in the Introduction. So, some spins placed far apart may interact, 
while pairs of spins placed relatively close to each might not. This can be of interest to 
model special disordered systems, having interactions of all ranges. 

For the packing problem, it is important to calculate the radius of each circle c(i, n), which 
depends both on the generation n as well as on the local environment, i.e., the generations of 
the three circles that it touches n p (i, k), k = 1, 2, 3. It is possible to include this dependence 
into the magnetic model, by defining coupling constants (or local fields) depending on the 
distance between nodes, i.e., J(n,i,n p (i, j)). We restrict ourselves, however, to a simpler 
situation, where the coupling constants J only depend on the generation n at which the 
edge was introduced into the network. In Figure lb we illustrate the first three steps of the 
construction of the model. 

To have a physically interesting model, it may be reasonable to choose values for J that 
increase with n. Indeed, when n increases, the average length of the edges introduced in 
that generation decreases and, as the spins get closer, we might expect that the interactions 
among them become stronger. However, to avoid the divergence of the energy within the 
lattice, it is necessary to renormalize the value of all J's as n increases. To accomplish this, 
we define J n ,m-, m = 1, ...,n as the value of the constant introduced in the m — th generation, 



when the lattice has been built up to its n-th generation, and require that J nim decreases 
with n — m. A sufficiently general choice would be 

(n — m + l) a 

where Jo may have ferro- (> 0) or antiferromagnetic (< 0) character, the exponent a controls 
how the interactions decay with the difference n — m, and yields a possibility of choosing 
the interactions according to the generation at which they were introduced. 

One of the extreme situations, a = 0, corresponds to equal interactions along all edges in 
the network. On the other hand, in the a — > oo limit, the model contains only finite inter- 
actions for the subset of edges that were introduced in that last generation n, as illustrated 
in Figure lb. The number of these surviving bonds is given by B(n) — B(n — 1) = 3™ _1 , 
so that the average coordination number is reduced to four. Moreover, the lattice is then 
composed by four-sided polygons, so that competition and frustration due to the presence 
of antiferromagnetic bonds can never occur in this a limit. 

To close this section we write down the formal Hamiltonian of this model 

H n = - JjjSjSj -hy^j, (2) 

where all pairs of nearest neighbors denoted by (i, j) are defined according to the construction 
rules of the network, and the constants Jjj must be chosen from the set defined in eq.(£Q) 
according to the value of m in which the edge was introduced. We also include a constant 
magnetic field h, which allows for the evaluation of magnetic properties. The notation in 
eq.(|2*|) does not include the selection of the bonds that are taken into account. The evaluation 
of a partition function can in fact be set up in much proper terms through the transfer matrix 
formalism. 



III. TRANSFER MATRIX AND RECURRENCE RELATIONS 

The numerical evaluation of the partition function for magnetic models on scale invariant 
graphs with a finite number of end nodes has been performed with the help of TM derived 
maps for a large number of lattices and models. The problem we are interested in this work 
is also suitable to be analyzed within this framework. For the sake of simplicity, let us first 
consider the homogeneous case a = 0, and set h = 0. If we consider the first generation 
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n = 1, we observe that a 2 x 2 TM Mi, which takes into account all interactions between 
the sites i and k of Figure 1 can be written as: 

la x b\ (a(a 2 + b 2 ) 2ab 2 \ 
Mi=\ = h (3) 

\b x a x J \ 2ab 2 a{a 2 + b 2 )J 

where a = b^ 1 = exp(/5Jo). Mi can be used to describe one single cell or a linear chain of 
triangles that are connected by their bases. It is also possible to define a 2 x 4 TM L\, that 
describes the interactions among sites i,j and k, where the column labels k are composed 
from the pair (j, k) according to the lexicographic order, i.e., n = 2(j — 1) + k according to 

(Ci d^ di c/i \ / a 3 ab 2 ab 2 ab 2 \ 
\=\ • (4) 

di di di c\ J \ab 2 ab 2 ab 2 a 3 / 

Of course we note that ai = c\ + d\ and b\ = 2d\. 

Within the proposed framework, all interactions between the sites i and k, for any higher 
order generations n = 2,3,4..., should be written in terms of a single 2x2 TM's M n , 
with the same distribution of matrix elements as Mi. Moreover, the matrix elements of 
M n should be written in terms of those of the matrices of the lower generation n — 1 only. 
This turns out to be feasible since the Apollonian lattice, in a generation n + 1, can be 
decomposed into three sublattices, each one of them being a deformed lattice of generation 
n. Since the coupling constants do not depend on the actual distance between the sites, 
each of the three sublattices entails the same coupling constants and magnetic structure as 
the n-lattice. Thus, a matrix M n+ i can indeed be written in terms of three matrices M n . 
To achieve this we remark that, in any generation, the three sublattices share their three 
outmost sites, which we label as i,j, k and £. This last one occupies the geometrical center 
of the n + 1-lattice. Then, M n+ i can be defined using 

(M n+ i)i t k = {Ln)i,j£{L n )j t £k (ify kjt- (5) 
U 

and 

{L n +l)i,jk = / J(-^n)i,^(-^n)i,tt(-^n)fc,3^- (6) 
i 

As one can easily observe by direct evaluation of eq.© and eq.(jOJ), all matrices M n and L n , 
share the same matrix element distribution as Mi and L\. So, it is possible to immediately 
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write down recurrence relations for the elements of L n+ i in terms of those of L n as 

C n +1 = Cn + dn 

dn-\-i c n d n -\- d n 

from which the elements a n+ \ = c n+ i + d n+ i and b n+ \ = 2d n+ i can be obtained. 

However, a direct evaluation of the matrix elements defined by the equations ©-(01 
shows that they do not exactly describe the interactions between the sites % and k. For 
instance, in the generation n = 2, the number of magnetic bonds is equal to the number 
of edges B{n = 2) = 6. On the other hand, we see that the Boltzmann weights in M 2 are 
expressed by combinations of exp(/5r J ) exp(— (3sJ ), with r + s = 9 instead of r + s = 6. 
This is due to the fact that each one of the interactions between the site I and its neighbors 
i,j and k appears twice in eq.©. To describe the thermodynamics of the system with the 
help of equations ©-(0, it is necessary to carry out a small correction, namely to redefine 
a and b as a = 6 _1 = exp(/3Jo/2). With this modification, the Boltzmann weights in each 
element are expressed by exponentials of /3(r — s) Jo, where r + s = B(n) — 3/2. So, for each 
spin configuration, the ratio between the correct energy and that one provided by equations 
© and © is roughly proportional to (-B(n) — 3/2)/ B(n), which — > 1 in the limit n — > oo. 

These definitions are sufficient for only the uniform interaction model a = 0. Further 
modifications in eqs.©-© are required to obtain the correct maps for general a. To cast 
this into a single recurrence relation, we first note that it is not necessary to use two labels 
n and m to insert the correct coupling constants into M n . As these matrices are recursively 
defined, the largest and most abundant J n , n corresponds to the constant introduced into 
M\, which is reproduced in ever growing number by the successive use of equations like 
©-©• On the other hand, the smallest and least frequent J n> x represents the constant that 
is inserted into the sequence of TM's exactly at the n-th generation. So we consider 

(±l) n J 

Jn,m — > Jn = 7^ _ ^ Q • (8) 

Note that the changes carried into the denominator of eq.© set Ji = 0, and J n +i — * J n - 
This strategy is necessary to avoid taking into account more than once the effect of bonds 
introduced when n = 1, and it is somehow equivalent to the redefinition of a and b discussed 
above. 
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Then we modify equations (J3J) and (0) according to 

{M n+ i)^k = '^^(L n ) i ji(L n )i^k(L t n )k,ji(C n )i/( y C n )ij(C n )e^ (9) 

and 

(L n +l)i,k — ^^(^n)i^(^n)i/fe(^^)fcj^(^n)i^(C'n)£,j(C , n)^,fc, (10) 

I 

where the 2x2 TMs C n are defined by 



c n = | 


' Vn 9n' 






H 









(11) 



exp(/3J n ) exp(-/3J n ) 
^exp(-/3J n ) exp(/3J n ) 

With these definitions, it is possible to observe that the number of non-zero coupling con- 
stants in the lattice is Bin) — 3 so that, in the n — > oo limit, equations (|5j)- (fTT|) accurately 
describe the thermodynamic properties of the model. The recurrence maps for the matrix 
elements derived from (fHJ)- (fTUj) read 

C n +1 = (? n Vn + ^n^n ^) 
d n +l — CndnPnQn + '^nPnQ'n 

From equations (7) or (11) it is possible to derive recurrence maps for the free energy 
/„ = — Tln(c n + d n )/N{n) and correlation length £ n = l/ln((c„ + d n )/(c n — d n )) at two 
subsequent generations, f n+1 = f n+1 (f n ,£ n ;T) and £ n+i = £ n +i(fn, £,n, T) can easily be 
derived |lRL llflj . This set of maps can be increased by working out explicit recurrence 
relations for the derivatives of f n (T) with respect to both the temperature and the magnetic 
field, obtaining the entropy s(T), the specific heat c(T), the spontaneous magnetization 
m{T) = m(T,h = 0), and the magnetic susceptibility x(T,h = 0). For this last purpose, 
we have to consider h ^ and insert it into the matrices M n . This modification breaks the 
up-down symmetry of the problem, so that the matrices M n and L n have a larger number 
of distinct matrix elements. This is a straightforward procedure that has been carried out 



for other models 
in this work. 



18l I19I ] . In the Appendix we present the full set of recurrence maps used 



IV. RESULTS 

We study the thermodynamic functions, i.e. the free energy /, the entropy s, the specific 
heat c, the spontaneous magnetization m and the correlation length £ as function of the 
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Figure 2 



FIG. 2: Thermodynamic functions for the ferromagnetic model. Solid, dashed and dotted lines 
indicate a = 0, 1 and oo. 

temperature T, as shown in Figures 2 to 5. They were obtained by numerically iterating 
the set of maps shown in the Appendix, starting with T dependent initial conditions, until 
a value n*(T). This temperature dependent value is set automatically, by requiring that 
one (or a set) of the intensive quantities and/or the correlation length, have converged to 
a fixed value, within a previously established relative tolerance. This is usually ~ 10~ 15 , 
as we work with double precision variables. Convergence based only on the value of / and 
its derivatives is much faster than for £, specially when the system is in the ordered phase. 
Otherwise stated, if we call n*(T),x = f,s,c,m,£, the value of n at which the function x 
has converged for that particular value of T, then we find that n|(T) always assumes the 
largest value. 

In Figure 2 a,b,c and d we show the entropy s, the specific heat c, spontaneous magne- 
tization m(h = 0) and correlation length £, for three distinct values of a, when all coupling 
constants have ferromagnetic character, i.e., J > and (+1)™ in eq.©. The qualitative 
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FIG. 3: Size dependent critical temperature for distinct values of a. 

behavior does not depend on the values of a(0, 1, and oo), i.e., whether interactions are 
only short (a — ► oo) or long range (a = 0). For all cases we see that, for low values of 
T, long range correlation sets in, as is evident from the spontaneous magnetization and the 
numerical divergence of £. The remarkable feature, however, points to the absence of any 
criticality when T is increased. The insert in Fig 2b shows that, when a = 0, m goes to zero 
smoothly, as exp(— T), with no evidence of a sharp transition to m = at a well defined 
critical temperature. If we consider a > 0, we still find a smooth, but stronger, decay, 
namely as m ~ exp(— T A ). The curve for the specific heat is also smooth, showing a typical 
Schottky maximum, again without any evidence of a divergence, that would be expected for 
a usual phase transition. 

The results for the correlation length (d) are also distinct from those found for other 
scale invariant models, as the diamond hierarchical lattice (DHL) 19]. There, £ is finite 
for large values T and numerically diverges for all values of T below a well defined critical 
value T c which in our case means attaining a value larger than 10 16 , the largest allowed 
number in our algorithm. Within this region, the actual value reached by £(T), has no 
precise meaning. Typically it is much higher than those in the disordered phase, and is also 
characterized by the presence of random fluctuations. As mentioned before, n^(T) is larger 
than n*(T),x = f, s, c, m, but even if we stop the iterations at n*JT), £ has already reached 
this very high plateau. This shows that it is not actually necessary to proceed further with 
the iteration of the maps, as we would obtain only a meaningless value for £. 

In the present case, if we use rij(T) to stop the iteration of the maps, we observe that £ 
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FIG. 4: Thermodynamic functions for the anti- ferromagnetic model. Solid, dash and dots indicate 
a = 0, 1 and oo. £ is drawn only for a = and reciprocal temperature axis. 

diverges at low temperatures, expressing long range order. When T is increased beyond a 
given value of T*, it converges to a well defined value, suggesting the break of long range 
correlation. However, if the iteration procedure is pursued to a value of n > n*JT) we 
observe that the T interval in which £ diverges becomes larger. This finding has driven us 
to proceed with the iteration of the maps in a different way. We fix a value n > n%(T = 1), 
and iterate the maps until reaching n for all values within a large T interval, as shown in 
Figure 2c. Then it is possible to precisely evaluate a critical value T c (n), as the value of T 
where the behavior of £ changes. In Figure 3 we show how T c (n) depends on n, for several 
distinct values of a. Our findings for this unusual kind of critical behavior suggest a power 
law T c (n) ~ n T ( a \ with r going continuously from r(a = 0) = 1 to r(a = oo) = 1/2 . We 
recall that a similar behavior has been reported for spin models on another scale-free lattice 

121- 

In Figure 4 we show that antiferromagnetic interactions (Jo < 0) change the thermody- 
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namic behavior of the model. The most interesting situation is observed for a = 0. All trian- 
gles in the lattice are frustrated and, as expected, a residual entropy so = 0.222 is measured. 
This best numerical value is smaller than that for the triangular lattice s = 0.3238.., jsjj], 
and much smaller than that obtained for the Ising model on the Sierpiski gasket Sq = 0.493.., 



At the same time we find finite well defined values for £ for all values of T, which are 
robust with respect to the value of n where the iterations are stopped. This is illustrated 
in Figure 4d, which also shows that, as T — > 0, £ decays like exp(— 1/T), typical for the 1-d 
chain. This is a somewhat unexpected behavior, as the presence of frustration usually does 
not allow for long range correlation of spin orientation, even at T = (e.g. the AF Ising 
model on the SG At the same time, this result must be related to the non- vanishing 

behavior for the magnetization, shown in Figure 4c. It indicates that the number of spins 
pointing in each direction is not the same. Once again, this behavior is different from that 
obtained for other frustrated lattices, like the planar triangular lattice or the SG. Finally, 
the behavior of the specific heat looks like those found when Jo > 0, for any value of a. 

In Figure 4 we also draw curves for the thermodynamic functions when a = 1, and oo. In 
Figure 4c we see that the magnetization curve always saturates at m(T — ^ 0, /z — > 0) = 1/3. 
This indicates that, in this limit, the number of spins pointing in opposite directions is 
not the same, as is the case for the triangular lattice, but stay in proportion 2/3 to 1/3 
independent on a. For some range of values of a, m goes through a maximum (e.g. at 
(T, m) = (1.2,0.351)) when a = 1, so that a reentrant behavior at low temperatures is 
observed. 

Let us also discuss how the presence of interactions with different sign affects the behavior 
of the system, i.e., when we take (— l) n in eq.(jHJ)- As expected, the result depends whether 
Jq > or < 0. In the first case, competition and frustration give rise to residual entropy 
when a = 0, as illustrated in Figure 5a. However we note that the value of so is smaller 
(~ 0.152) than in the case of equal AF interactions. This happens because not all triangular 
units are build by an odd number of AF bonds, as we can easily see by inspecting the first 
generations with the help of Figure lb. Also for a = we note the presence of a double 
Schottky peak in the specific heat (Figure 5b). For q^0 there is no remarkable difference 
between the curves for sore with respect to those obtained for interactions with the same 
sign. The results also show that the low temperature magnetization saturates at the value 
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FIG. 5: Thermodynamic functions for the model with alternating ferro and anti-ferromagnetic 
coupling, starting with Jo = 1. Solid, dash and dots indicate a = 0, 1 and oo and Jo > 0. In (c) 
three curves for the magnetization when Jo = — 1 and same values of a saturate at m = 1/3 are 
also drawn. 

m = 7/9 only when a = 0, otherwise m — 1. 

Finally, when J < and alternating sign are considered, no frustrated bonds and, 
consequently, no residual entropy is found. The curves for the specific heat are also smooth 
like all other cases. The magnetization curves saturate again to m = 1/3 as T — > for all 
values of a. However, reentrant behavior similar to that found for some of the AF cases has 
not been observed, so that the typical shape is that shown in Figure 5c. 



V. CONCLUSIONS 



We have studied a family of Ising models on an Apollonian network using the transfer 
matrix technique. On one hand we considered ferro- and antiferromagnetic couplings and 
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on the other hand we generalized the interaction as being dependent on the generation as 
well in the sign as in the strength quantified by an additional parameter a. 

For purely ferromagnetic couplings we always find order in the thermodynamic limit 
independent on a which is in agreement with what has been found on other scale-free 
lattices^?}]. Interestingly the effective critical temperature at which the correlation length 
diverges goes to infinity with the system size with a power-law in the number of generations 
with an exponent that depends on a. For antiferromagnetic couplings we find a disordered 
phase for any finite temperature but a diverging correlation length at T = 0. This later 
observation is unusual as it does not appear for instance on the Sierpinski gasket j^. 

Considering the Apollonian network as a model for the connections between cities as 
described in ref.[14j our result can be applied to the formation of opinions where spin up 
means one opinion and spin down the other one. The results for the ferromagnetic case 
implies that independent on the strength of the couplings between the cities as long as it is 
not zero one single opinion will finally prevail. 

If the Appolonian network describes the force lines in a dense polydisperse packing with 
each particle having a magnetic moment as it is the case in tectonic faults |24j our result 
for the ferromagnetic Ising model would imply that if all particles have a moment of equal 
strength one always finds a spontaneous magnetization. 

Our calculations can be generalized to random couplings (spin glass) which in fact is 
work in preparation. One can also imagine studying other more complex magnetic models 
on the Apollonian networks, like the Potts model, the XY model or the Heisenberg model 
and one can also study the magnetic properties of Apollonian packings of different topology 



or higher dimension and even the case of the random Apollonian packing 



25|. 
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VII. APPENDIX 



The maps for the free energy and correlation length derived from eqs. (|T2"j) read: 



fn+i = - ]v^{3 lna n + ln{l + 3z n (3 n (2 + (3 n ) + 3^(1 + 2&(1 + /3„) 2 ) + 



! + 

z3(l + 2/3 n )(2 + /3„) 2 }-61n2}, 



(13) 



(14) 



Cn+l — £n{l + £ra 

{ln(l + 3z n /3 n (2 + /3 n ) + 3^(1 + 2/ 3 n (l + f3 n f) + ^(1 + 2/3 n )(2 + /?„) 2 )- 
M*- 1 /^ + ^(1 + 4/3 n + 2^ + 2(31) + 4(4 + 10P n + 11/3= + 2/3*)+ 
4(3 + 10/3 n + 10/3,2 + 4/53))}}- 1 
where 2 n = (c n -d n )/(c n + d n ), a n = p n + q n and (3 n = (p n -q„)/(p n + q n ) and N n = N(n 
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